#include "lombscarglecalculator.h"
#include <cmath>

LombScargleCalculator::LombScargleCalculator(QObject *parent) : QObject(parent)
{
    // 构造函数
}

// 计算 Lomb-Scargle 周期图
QPair<QVector<double>, QVector<double>> LombScargleCalculator::Lombcalculate(const Eigen::VectorXd& t, const Eigen::VectorXd& y, int nfft, double maxF)
{
    QVector<double> PSD(nfft);
    QVector<double> F(nfft);

    double deltaF = maxF / nfft;
    for (int i = 0; i < nfft; ++i) {
        F[i] = i * deltaF;
    }
    PSD=calculate(t,y,Eigen::VectorXd::Map(F.data(),  F.size()));
    return qMakePair(PSD, F);
}

// 计算 Lomb-Scargle 周期图
QVector<double> LombScargleCalculator::calculate(const Eigen::VectorXd& t, const Eigen::VectorXd& y, const Eigen::VectorXd& f)
{
    int N = y.size();
    Eigen::VectorXd z = y.array()  - y.mean();

    QVector<double> Pn(f.size());
    for (int i = 0; i < f.size();  ++i) {
        double w = 2 * M_PI * f(i);
        if (w > 0) {
            Eigen::VectorXd twt = 2 * w * t;
            double tau = atan2((twt.array().sin()).sum(),  (twt.array().cos()).sum())  / (2 * w);
            Eigen::VectorXd wtmt = w * (t.array()  - tau);

            double cosTerm = (z.array()  * wtmt.array().cos()).sum();
            double sinTerm = (z.array()  * wtmt.array().sin()).sum();
            double cosSq = (wtmt.array().cos().square()).sum();
            double sinSq = (wtmt.array().sin().square()).sum();

            Pn[i] = (cosTerm * cosTerm / cosSq) + (sinTerm * sinTerm / sinSq);
        } else {
            Pn[i] = (z.array()  * t.array()).sum()  * (z.array()  * t.array()).sum()  / (t.array().square()).sum();
        }
    }

    for (int i = 0; i < Pn.size();  ++i) {
        Pn[i] /= N;
    }

    return Pn;
}
